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ABSTRACT 



This thesis investigates the pressure drag coefficient in 
the transonic regime over an axi- symmetric body, with a set of 
unique contour surfaces developed in a previous thesis. The 
contour surfaces were obtained by an exact solution of the 
small perturbation transonic equation, using the guidelines 
and tools developed at NPS . In this work. Computational Fluid 
Dynamics (CFD) was not only used to compute the afterbody 
contour surface, but also to investigate a conical afterbody 
and complete bodies, which are composed of an arbitrary 
forebody (ellipsoid) and variable afterbody (contour and 
conical) . Euler as well as Navier- Stokes flow- solvers were 
applied to the geometries of interest, giving Mach-number 
contours for viscous and inviscid flow, pressure drag 
coefficient magnitude, and depicting shock wave locations. On 
the basis of these results, it can be verified that our 
contour surface afterbodies will decrease by 15% the peak of 
the pressure drag coefficient (C^) versus Mach number curves 
in the transonic regime. These results can be used to design 
low pressure drag surfaces for such shapes as missiles, 
projectiles and aircraft engine nacelles. 
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I . INTRODUCTION 



The flight of an object over a wide range of speeds has a 
critical transition zone where both subsonic and supersonic 
types of flow exist. This speed regime is referred to as the 
transonic range. The critical aerodynamic behavior occurs in 
the range 0.8 < M < 1.2 depending on the object, where the 
aerodynamic coefficients have been found to change by as much 
as 100%. The behavior of the aerodynamic force components is 
usually characterized by a rapid increase in the coefficients 
followed by a sharp drop, in other words, a peak value for the 
coefficients arises in this regime [Ref.l and 2]. 

One of the aerodynamic components, aerodynamic drag, 
represents a significant adverse force on all flying objects 
such as aircraft, missiles and projectiles. A high drag force 
reduces the craft's range capability or equivalently requires 
more energy to achieve a certain range. Any effort to reduce 
the drag coefficient in the design process must concentrate on 
reducing the wake and pressure drag (inclusive of wave drag) 
contributions to the total drag. 

A traditional approach to investigate aerodynamic 
characteristics is based on wind tunnel data and actual flight 
testing as a source for improvement of configurations to get 
the best results. Unfortunately, both wind tunnel and flight 
tests are considerably expensive and time consuming. 
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particularly in the transonic flow regime. Because of the 
great lateral influence of profiles at transonic speeds (shock 
reflection of the wall bounces back to the model) , models must 
be made extremely small compared with wind tunnel dimensions 
and this introduces great experimental difficulties. 
Consequently, the data obtained at transonic speeds is 
considerably less reliable than at either subsonic or 
supersonic speeds. In addition, from a mathematical point of 
view, even the two dimensional small perturbation potential 
equation for transonic flow retains one non-linear term which 
is essential for non-divergent solutions at Mach one, but this 
non-linear partial differential equation has proven to be 
difficult to solve. Such inherent difficulties, coupled with 
the presence of shocks in the flow which cause boundary layer 
separation, have resulted in the creation of many approximate 
methods of solution which are employed in the design of 
transonic airfoils and the like [Ref .3, 4 and 5] . On the other 
hand, the use of numerical simulation known as Computational 
Fluid Dynamics (CFD) to predict aerodynamic characteristics 
greatly increases possibilities to improve design optimization 
at relatively low cost and allows for ease of design changes. 

Finally, using the latest capabilities of Euler as well as 
Navier- Stokes flow- solvers , it has been possible to compute 
the flow over axi- symmetric bodies with various contours in 
the transonic regime. 
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II TRANSONIC FLOW 



Transonic flows are characterized by the simultaneous 
presence within the flow field of both subsonic and supersonic 
regions. The properties of transonic flows can be derived from 
the general equations of gas dynamics, namely, the equations 
of state, continuity, momentum and energy. The following 
derivation based on assumptions that the flow is steady, 
irrotational and isentropic with no energy transfer, no body 
forces and no shear stresses (inviscid flow) . 



1. Small Perturbation Theory 

Starting with the equation of motion for steady, 
isentropic, inviscid flow in the index notation form [Ref .4] : 



du^ du^. 

ax, ■ ^ dx. 



( 1 ) 



Small - Perturbation Theory gives the velocity field as 



= U + u ; = V ; = w (2) 

where U denotes the free stream velocity in the x direction 
and u, V, w are called perturbation velocity components in the 
direction x, y, z respectively. 
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Writing eq. (1) out in full, and substituting the velocity 
field eq. (2) , gives the equation in terms of perturbation 
velocities . 



,2/au ^ ^ = 



\ dx dy 



dz 



dx 



dv 

dy 



dw 

dz 



( 3 ) 



+ (U^u)vlp * * vJ|i: * |5f) * . |h 



\dy dx] \dz dy 



\ dx dz 



From the energy equation for a perfect gas, the speed of sound 
(a) can be expressed in terms of the perturbation velocities. 



( U+u) 2 



(Y-1) 



_£2 

2 



(y-1) 



( 4 ) 



or 



,2 = 



(Y~1) 

2 



(2uU + + w^) 



( 5 ) 



Substituting eq.(5) into eq.(3), dividing by a„^, and 
rearranging the terms, gives the full exact equation in terms 
of perturbation velocities and free stream Mach number. This 
equation contains linear terms on the left-hand side, but on 
the right-hand side the terms are nonlinear. 
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(1-Mi) ^ ^ ^ ^ ^ = 
ox oy oz 



m: 



iy + 1) 



+ m: 



+ m: 





(y+1) 


i±- + 


(y- 


1) 


( v^+w^) 


du 


u 


2 


f/2 


2 




f/2 


dx 


R + 


(Y + 1) 


+ 


(Y- 


1) 


(w^ + u^) 


dv 


u 


2 




2 






dy 


R + 


iy*i) 


r1 + 


(y- 


1) 


{u^ + v^) 


dw 


U 


2 




2 






dz 



37 . ( 6 ) 



+ Ml 





fin. 


_dv\ 


+ ^/i + “\/ + 




+ RRl 


dw 




[w u)i 


[dy 


dx) 


WMd^ 


dx) 


uA 


Jy 


dz) 



If the perturbation velocities are small (u/U,v/U,w/U << 1) , 
eq. (6) can be simplified by neglecting the terms containing 
squares of the perturbation velocities on the right hand side, 
yielding : 



\ -*W / -N ^ 

OX oy oz 



U dx U\dy dz 






U\ dy dx 



U\ dz dx 



( 7 ) 
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For further simplification, in eq. (7) all the terms on the 
right hand side can be neglected, in comparison to those on 
the left hand side. This gives the linear equation, which 
contains only perturbation velocities and is valid only for 
subsonic and supersonic flow. 



(1-Mi)#; 

ox 



dv ^ 
dy dz 



= 0 



( 8 ) 



For transonic flow, where -> 1, the coefficient of du/dx on 
the left hand side becomes very small, but it is not correct 
to neglect the first term on the right hand side of eq. (7) . 
Therefore, the governing equation of transonic flow in term of 
perturbation velocities is as follows: 






du 



dv 



dw 

Ti 






( 9 ) 



For irrotational flow, a perturbation velocity potential 0 
exists , 



u 






( 10 ) 



Substitution of eq. (10) into eq. (9) , gives the governing 
equation for transonic flow in terms of the velocity 
potential . 
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( 11 ) 



(7-^1) dip d^<p 
dx^ dy^ dz^ ^ dx^ 

For bodies of revolution, it is convenient to use cylindrical 
coordinates (x,r,0) where x is aligned with the body axis. The 
velocity components corresponding to (x,r, 0 ) are u,^, and 
U 0 , respectively. The velocity potential in the cylindrical 
coordinates , 

- TT - ^0 _ d(b _ 1 dd) /io\ 

Transforming Cartesian coordinate eq.(9) into cylindrical 
coordinate, gives the governing equation of transonic flow. 



(1-M^) ^ ^ ^ d<p d^4> (13) 

dx^ dr^ r"5r qq2 u "5x 

For axially symmetric flow, where the conditions are the same 
in every meridian plane, there is no variation with 6, so the 
small perturbation, non-linear, axi- symmetric transonic 
potential equation can be written as follows ; 



d^d> ^ 1 d<p ^ H 

dr2 r~5r Qyr 2 u "53c ^^2 



(14) 
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Or rewriting in shorthand notation, 




( 15 ) 



Reference 6 introduces a modified potential equation for axi- 
syitimetric flows 



2. An exact solution for axi- symmetric, transonic flow. 

Solutions to the modified transonic eq. (16) have been 
given by Biblarz [Ref. 6 and 7] by using the separation of 
variables approach with a potential function 4>{x,r) of the 
form 



Substituting the above 0 function in the modified transonic 
eq. (16) results in two ordinary, second order, non-linear 
differential equations 




(16) 



where the modified velocity potentials are : 



^ wf (7+1) ^ (7+1) 



(17) 



4>{x,r) =|(x)7?(r) + (l-Mf)x 



(18) 
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(19) 



dx dx2 



- X? = 0 



and 



d^r\ 

~d^ 



1 dt) 
r dr 



- Xr?2 = 0 



where X is Che separation constant . 

The solution to the first differential eq. (19) is 
by multiplying both sides d^/dx, 



or 



^ d| d^g 
dx dx (^2 




(X|) = 0 



A HA? - 

dx 3 \ dx/ 2 ^ 



= 0 



Thus 



1 




where a is a constant given by Ref. 6. 
Rearranging eq.(23) 



dx = 




Q! 



1 

1 



( 20 ) 



obtained 



( 21 ) 



( 22 ) 



(23) 



(24) 
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and integrating eq. (24) gives. 



X - Xq 




( 25 ) 



where a and Xq are integration constants. 

The solution to the second non-linear ordinary 
differential equation (20) , is obtained by an outer expansion 
method. [Ref. 8] 

Vir) = + (1-M,^)fi(r) + (1-M^)2f2(r) + ... (26) 

Where (l-M^^^) represents a small parameter and the first term 
is the purely sonic solution. 

By taking first and second derivatives eq. (26) and 
subtituting them into eq. (20) yields, 



1 ? (r) =. 



Xr^ 



/ 1 \ 2 

4- I 1 -W^ I aXr ^Ii_^a 2 x 2^2 (v^F. 2 ) 

no o*/o 



(27) 



Al-Hashel [Ref. 9] reported on eq. (26) and eq. (27) , 
implemented the boundary conditions with the constants a, 
a and X, gave the final results in the new variables. 



f - I 



3X 



2C-, 



(28) 
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_f(/8-+2) = ^(yS-+2) ^ 



( 29 ) 



X = C-, 



1 1 
6 J 3\\ 2 






(30) 



1? = 



1? X 



(a X) 



0.4142 



(31) 



then 



X = 



= I 



d| 






l-M, 



2 I 1.7574 



(32) 



where Xq=0, and 







£=(v/F+2) 






50 . 63 



(s/'B’+2) 



(33) 



Equation (32) has been numerically integrated and plotted on 
figure (1) as ^ (x) versus x, and eq. (33) will be evaluated and 
plotted in figure (2) as ^ (f) versus f for M^, = 1.05, 1.1 and 
1.2. A "patching" technique discussed in Ref. 9 has been used 
here . 
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III. PRESSURE COEFFICIENT AND BOUNDARY SURFACES 



1. Pressure Coefficient 

Liepmann and Roshko [Ref .4] define the pressure coefficient 

as , 



(P-Pco) 




( 34 ) 



From the isentropic relation, we have pressure ratio in terms 
of Mach number and after substituting into eq. (34) yields 





O 


Y/(Y-1) 


2 


2 + (7-1)AC 


- 1 




2+(y-l)M^ 





( 35 ) 



Introducing = U^/a„^, = u^/a^ and using energy equation 
(recall eq.4 and eq.5), the pressure coefficient can be 
expressed in terms of the perturbation velocities, 



C„ = 










r ^ 




y/(y-i) 


« 


L+ ”^ ^ 


^ _ iU+u) 2+v^+w^ 




> 

- 1 




2 “ 









( 36 ) 



2 

yM^ 



2 “ 



2u u^+v^+w^ 



Ty/(y-i) 



U 






1 
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Using the binomial expansion on the expression inside the 
square brackets in eq. (36) , we obtain the pressure coefficient 
in the form 



= - 



. v^ + w^ 



u c/2 



( 37 ) 



For axi- symmetric flow, in cylindrical coordinates where u=u^ 
and (v^+w^) = u,.^ , substituting into eq. (37) yields 



= - 



2 Liy. 2 / 


^xY 


UrY 


+ 1-Afi 


— ^ 








J 



( 38 ) 



The linearized pressure coefficient approximation for axi- 
symmetric flow turns out to be 



C = - ^ 

c/„ 



( 39 ) 



Recall the modified axial velocity potential, eq. (17) , 



u. 



4 >^ = m ;(7 + 1 ) 



( 40 ) 



thus 



^ ^ 4>x 

mJ(7 + 1) 



( 41 ) 
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Substitute eq. (41) into eq. (39) , yields 



, _ ^ 

P — 2 



( 42 ) 



The derivative of the potential function eq. (18) with respect 
to X becomes 

+ (1-wf) (43) 

Ref . 6 introduces relation of the constants oi, a and X in 

expression. 



a = ±Ci|l-Mf ! 



1.7574 



(44) 



and 



'Cil ^1.242e 
"X" xO.7574 



1. 08x2 0 



(45) 



Rewriting eq. (23) and substitute the constant a eq. (44) 
becomes 



dx 



3\ v2 



+ Cl I (l-M^) I 



1.7574 



1 

3 



(46) 
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Recall Che new variable eq. (28) 



I = e 



3X 



2C-, 



(47) 



Substituting eq. (47) into eq. (46) and factoring out, yields 



1 




(1-m:) 



1.7574 



(48) 



Recall and arrange new variable eq. (31) 



V 



(a X)0.4142 

X 



7 



(49) 



Rewriting eq. (43) Che modified potential function and 
inserting eq. (48) and eq. (49) , becomes 



K - » 



Cl [f- 



(1-Aff) 



1.7574 



(l-M^) 



(50) 



or 



4>x = 



1 

^ 0 . 4142^.3 

V 



X 



0.5858 






(51) 
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Substituting the expression of constants eq. (45) into eq. (51) 



= 0.2208 rj p + \ (1-M^) 11-7574 



(1-wf ) 



( 52 ) 



Finally, rewriting the pressure coefficient as 



= 



-2 



(7 + 1 ) 



O.2208^[p- 



ii-Mh 



1.7574 



^ + (i-wf ) 



( 53 ) 



2 . Boundary Surfaces 

For an inviscid flow, the condition to be applied at the 
surface of a solid boundary is that the direction of the flow 
velocity vector must be tangent to the solid surface. In other 
words, the velocity vector is everywhere at right angles to 
the normal to the solid boundary [Ref. 10] . In addition the 
boundary condition requires that the gradient of potential p 
vanish far ahead of the body. In terms of perturbation 
velocities this boundary condition becomes 



dx 



suiface 






( 54 ) 
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Recall the modified velocity perturbation potential eq. (17) 



u. 



4>r = ( 7 ^ 1 ) 






(55) 



thus 



Ur ^ <t>r 

wf (7+1) 



(56) 



substituting eq. (56) into eq. (54) , 



l—\ = 

wf (7+1) 



(57) 



Taking the differential of eq. (18) with respect to r gives 






V ^ 

^ dr 



(58) 



Expressing eq. (58) in terms of new variables, recall eq. (28) , 
eq. (29) and eq. (31) 



(pr = 0.0848 ? 



drj 

'df 



(59) 
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Recalling eq. (33) and taking the derivative with respect to f , 
then substituting eq.(59) into eq. (57) , yields 





0. 


,0848 ^ 


— +2.8284 \l-M^\ f ^-®284+0 . 1512 (1-Wi) 2 f6.657 


\ dxl 




(y+1) 


^ II AW# 

. ^ 



( 60 ) 



Replacing the left hand side of eq. (60) with the new variables 
eq. (29) and eq(30), then rearranging it becomes 



dx 



-+2.8284 |l-wf| f ^-^284+0 . 1512 (1-Wf) 2 

Afi(Y+l) 



( 61 ) 



For further arrangement, by separating variables eq. (61) , 
yields 



df 

— -2.8284|l-Afi| f ^-®284 _o . 1512 (i-Afi) 2 
f ^ 



-0 . 0326 

MUy*i) 



I dx 



( 62 ) 



Then, integrating both sides of eq. (62) becomes 
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( 63 ) 



J 



— -2. 8284ll-Mi|f -0 . 1512 (1-wi) 2 _f6.657 

f ^ 



fd(Y*i) ■i 



where 



1.2074 

11 -^ 2 | 0.2071 



( 64 ) 



Al-hashel [Ref. 10] has developed and computed eq. (63) 
using numerical integration, to determine the boundary- 
surfaces in dimensional and non-dimensional (normalized) form 
for M„ = 1.05, 1.10 and 1.20 as depicted in Figure (3) and 
Figure (4) . Based on these calculations, the geometric grid of 
the afterbody for Mach number 1.10 and 1.20 are developed for 
further study with Computational Fluid Dynamic (CFD) . 

This thesis research also examines conical afterbodies as 
a solid afterbody boundary surface with base diameter ratio 
(djj/djj,) of 0.50 and 0.75 and conical turning angle (/J) of 26.6 
and 14.0 degree respectively [Ref. 11 and 12]. Then, for 
further investigation, the complete bodies as a solid boundary 
surface are generated, with a kind of forebody (ellipsoid) 
joined with the contour surface afterbodies as well as the 



conical afterbodies. 



IV. DRAG 



Drag is one of the aerodynamic force components parallel 
to free stream velocity (Uj^) . It represents a significant 
adverse force on all flying objects. Basically, the drag force 
is divided into two categories, the drag caused by forces 
acting normal to the boundary surface which is called pressure 
drag (inclusive forebody, base and wave drag) and that arising 
from the tangential forces acting on the surface, by virtue of 
viscosity which is called viscous drag or skin friction. 
Figure (5) shows the components of the drag as function of 
Mach number and the methods used to compute body drag in four 
Mach number regimes [Ref .5] . 




Figure 5. Components of the Drag [Ref. 5] 
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1. Skin Friction 



Skin friction is the result of shearing force within the 
boundary layer of a viscous flow, acting tangentially to a 
surface in motion relative to the fluid. The amount of viscous 
resistance depends upon whether the flow is laminar or 
turbulent [Ref. 5] . Krasnov [Ref .13] introduces the laminar and 
turbulent skin friction for flat plates, based on the boundary- 
layer theory, for compressible flow. These formulas also valid 
for bodies of revolution with infinite length. 

Laminar flow. Re < 10® : 



^fL 



(1 + 0.03Ad) ^ 

{Re) 2 



(65) 



Turbulent flow. Re > 10® : 



'fr 



[logio-Se 

(1+0. iwf) ^ ^ 



( 66 ) 



where Re is the Reynold number 
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2. Pressure Drag 

Liepmann and Roshko [Ref. 4] introduce the pressure drag 
formula (inclusive wave drag and base drag) for axi- symmetric 
bodies as, 



L 

D = f p dS - PgSiL) 

0 

(67) 

L 

= J (p-pJdS+ {p^-Ps)S{L} 

0 



or, in dimensionless form, 



D 

Q^SiL) 



1 

S(L) 



f C 

J Pdx 



dx + C, 



pB 



( 68 ) 



where S (x) = irr^ is the cross-sectional area of the body at x. 



dS{x) 



2Tirdr = 2izr-^dx 
dx 



(69) 



or 



dS 

dx 



2nr 



dx 



(70) 
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'' The first term in square bracket eq. (68) can be solved by 

I inserting eq.(53), eq. (60) and eq. (70) yielding, 

I 



L 


L 




r 





< 


-2 n 


P~d5t 

0 


t 


i 


Wf (7+1) ■- 






(71) 



(27rr) 



M„(7+1) 



0.03 26 1 1—8 ^2. 8284 | 1-mJ | r'^ +0 .1512 (1-wf ) 2^6-657 



f3 



dx 



Recalling new variables eq. (28) , eq. (29) , eq. (30) and 
inserting = 1-Mo,^ into eq. (71) yields, 



I 






ri [0.22087J [| + | 

i (1-/32) (7+1) 




1 

5. 13 02 Cl® 
j^O.7071 ^0.2071 



(72) 



0.0326 
(1-/32) 



•I 



-8 



+2 . 82 84 I /32 I • ®284 



|dx 



where C^, a and X are constants. 
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The second term in square bracket eq. (68) is the contribution 
to the drag by the base pressure pg. Ref. 4 claims, that the 
values of the base pressure coefficient Cpg must be obtained 
experimentally. However Krasnov [Ref. 13] introduces the 
boattail drag coefficient for a conical afterbody as. 



CpB = 0.002 



0.8 + 



Ml 






( 73 ) 



where : 

6 = turning angle of the conical body 

Sg = ratio of the base area to the mid- section area 

In this work, the skin friction coefficient as well as the 
pressure drag coefficient of the geometries of interest, at 
specified Mach number, can be obtained directly from the 
results of CFD calculations. The Euler solution will give the 
pressure drag and the Navier-Stokes solution will give both. 
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V. COMPUTATIONAL FLUID DYNAMICS 



The rapid advancement in the speed of computers and their 
enormous memory size has led to the emergence of the field of 
computational fluid dynamics (CFD) . This branch of fluid 
dynamics complements its experimental and theoretical 
branches, by providing an alternate, cost-effective means of 
simulating real flow. It also offers the only means of 
examining theoretical advances for conditions unavailable 
experimentally. As a model based method, CFD can provide the 
convenience of being able to switch off specific terms in the 
governing equations [Ref. 14], so as to assist the researcher 
in understanding the contributions of various physical 
factors . 

In this work, CFD was used to compute the axi- symmetric 
flow over the afterbody geometry of models only (the boundary 
surfaces obtained by the small perturbation method [Ref. 9] and 
conical afterbody) and over complete body models which are 
composed of forebody (ellipsoid) and afterbody. 

1. Grid Generation 

The computer programs GRAPE [Ref. 15] and GRIDGEN2D 
[Ref. 16] are tools used to generate two-dimensional structured 
grids about airfoils and other shapes by the use of algebraic 
or Poisson differential equation solvers. GRAPE was used for 
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the geometry with smooth contour surfaces, while GRIDGEN2D was 
used for the geometry with a conical afterbody with non- smooth 
points. Outer and inner boundaries were specified as the C- 
type grid for afterbody models only, while the 0-type grid was 
for the complete body models, where both type of grids treat 
the surface of the body as the inner boundary. The important 
characteristics in a grid generation technique are the ability 
to specify the spacing between mesh points at the boundary, in 
the direction normal to the boundary, and the control of the 
angles with which mesh lines intersect the boundaries which is 
known as orthogonality [Ref. 15]. 

Figure (6) is a typical output of program GRAPE for an 
afterbody with small perturbation solution contour for Mach 
1.10 (SPS_1.1), C-type grid with the grid size of 115x60, and 
figure (7) is the complete body with the grid size of 152x60 
(0-type grid) . 

To develop a three-dimensional grid from two-dimensional 
grid (output program GRAPE), the FORTRAN code called HALF.F 
(Appendix) was used to write out the half of the 2D grid and 
convert them into 3D plane -grid. Then, the FORTRAN code 
ROTATEGR.F (Appendix) was used to rotate the 3D plane-grid for 
11 planes plus 2 more to generate a 3D volume grid of half of 
an axi- symmetric body surface as shown in figure (8) and 
figure (9) with grid size of 58x13x60 and 77x13x60 
respectively. These grids are ready for further processing 
with a flow- solver. 



26 



Program GRIDGEN2D gives the pattern of how to build 20- 
grid. It consists of four subfaces which are treated as 
boundaries. Grid spacing is determined by setting up the 
distribution points of each pair subfaces. It uses equal 
spacing along the body surface, while the direction normal to 
the body surface uses geometric spacing with specified width 
in the beginning. Figure (10) and Figure (11) are typical 
outputs of GRIDGEN2D for an afterbody and for a complete body 
with the grid size of 58x60 and 77x60 respectively. FORTRAN 
code D2D3.F (Appendix) is used to convert 2D grid output from 
GRIDGEN2D into 3D plane -grid. Then we apply the FORTRAN code 
ROTATEGR.F to generate a 3D volume grid of half of an axi- 
symmetric body as depicted in Figure (12) and Figure (13) with 
the grid size of 58x13x60 and 77x13x60. 

In this research, we also attempted to develop a fine grid 
for complete bodies (0-type grid) , where the radius of outer 
boundary is set up to be five times of the body's length. The 
steps to obtain 3D volume grid are similar to those described 
above with the final grid size of 77x9x120 and the final 
geometry of a quarter of an axi- symmetric body instead of a 
half body. Figures (14) and (15) are typical fine grids of the 
complete body with a conical afterbody and a complete body 
with contour surface afterbody. 
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2 . 



Flow- solver 



The OVERFLOW program [Ref. 17] was developed by NASA Ames 
Research Center. It uses either 3-D Euler or Navier-Stok.es 
flow- solvers for inviscid/viscous flow, by setting True or 
False the parameter VISINP (viscosity input) in the input file 
(overflow. in) . Before flow-solver code (OVERFLOW) is applied 
on the grid file, a formatted 3D grid file named "grid. for" 
must be converted into an unformatted grid file named 
"grid. in", by using FORTRAN code called READX.F (Appendix). 
Then, the NAMELIST input file parameter specification must be 
written for running OVERFLOW, it is called "overflow. in" . The 
input parameter consists of the number of iterations, 
timesteps, calculation methods, smoothing, type of flow and 
boundary conditions for each grid. The value of angle of 
attack (ALPHA) depends on the orientation of the grid in the 
coordinate system. In this case, ALPHA is 180° (flow comes 
from the x-positive to the x-negative direction) . The boundary 
conditions depend on the geometry corresponding to the final 
3D volume grid. Both input parameters, the angle of attack and 
the boundary conditions for each geometric shape, are 
tabulated in Table l. 

The file overflow, in (Appendix) is a typical input 
parameter specification of the axi- symmetric body with a 
conical afterbody using the grid size of 77x13x60. This file 
input uses NSTEP=100, Mach number=0.95, ALPHA=180°. The 
calculation method depends on a central difference Euler term 
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in J, K, L and ARC3D diagonal factorization and dissipation 
scheme (IRHS = 0, ILHS=2 , IDISS=2) . For the first attempt a time 
step DT=0.1, ITIME=1 and CFLMIN=5 . 0 was utilized. The boundary 
condition consists of 6 boundaries (NBC=6) , with the type of 
BC IBTYP=30 (inflow) in the J positive (1) , IBTYP=15 (Axis K 
round) in J negative (-1) direction, IBTYP=12 (symmetry in Y) 
in the K positive (2) and K negative (-2) direction, IBTYP=1 
(Inviscid adiabatic wall) in the L positive (3) direction and 
IBTYP=32 (Supersonic /subsonic inflow/outflow) in the L 
negative (-3) direction. 

The OVERFLOW program gives output files such as ovr.out, 
q.save, resid.out and fomo.out. To verify that a calculation 
is appropriate or converged, it can be traced by looking at 
the plot of residual history (resid.out). A convergence 
criterion was defined as the reduction in residuals by two 
orders of magnitude. If the first run (NSTEP=100) does not 
fulfill the convergence criterion, one may perform a restart 
by further running OVERFLOW and changing the parameter 
RESTART=.T. in the input file (overflow . in) and by copying 
q.save into q. restart. Repeating the above steps until all 
convergence criteria are met. If in the output file q.bomb 
appears, one tries another run by changing parameter time step 
DT and CFLMIN until the criteria are met. Then, one converts 
the output file q.save (unformatted) into formatted file named 
q.form by using FORTRAN code READl.F (Appendix) . 

Finally, the last step is to run PL0T3D program with the 
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source formatted files x=grid.for and q=q.form. By using 
functions on the PL0T3D program, many plots such as pressure, 
velocity, Mach number, vorticity, etc., are obtained. 

3. Results 

The results of CFD programs are grouped into the 
corresponding geometric shape, namely, afterbody only and 
complete body. The Euler flow- solver (inviscid flow) was 
applied to all axi- symmetric bodies, except for the afterbody 
models only, where both Euler and Navier- Stokes (viscous flow) 
were applied. Most of the calculations converged in 500 
iterations, meaning that the residual history achieved a two 
order of magnitude drop . 

Results can be analyzed by plotting the q.form file, 
output file from the OVERFLOW program, using the proper 
function in PL0T3D program (Mach number) . By interpreting the 
Mach number contour surrounding the body surfaces , one can be 
determine the characteristics of the flow field. In addition, 
the drag coefficient C^, can be obtained in the file fomo.out, 
output from OVERFLOW. Hence, in the sequences of Mach number, 
one can describe the significant flow characteristic of each 
geometric shape. 

a . Af terbodi es 

An Euler as well as a Navier- Stokes flow- solver was 
applied to these afterbody models. The approaching free stream 
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Mach number (ranging from 1.05 to 1.50) starts from the mid- 
section of the afterbodies, then the flow follows along the 
afterbody surface until it reaches a maximum local Mach number 
as tabulated in Table 2 . The maximum local Mach number for 
viscous flow is always lower compared to in the inviscid flow. 
This may be caused by the viscous flow itself since we are 
taking into account the shear force in the boundary layer near 
the surface. 

For the afterbodies from the small perturbation solution 
contour (SPSl.O, SPSl.l and SPS1.2), in inviscid flow, the 
shocks are formed at the contour surface. The location of the 
shock depends on the specific afterbody contour and the 
approaching Mach number; at the higher Mach numbers the shock 
appears a bit further downstream as shown in Figures (16) 
through (21) . As shown in Figure (22) and Figure (23) , this is 
apparently a result of viscosity; it shows the boundary layer 
by the increment of Mach number away form the surface. Weak 
shocks were formed further downstream compared to inviscid 
flow. Flow separation occurs in the starting contour region 
and is followed by circulating flow in the base region. 

For the conical afterbody, the approaching free stream 
Mach number increases following the mid- section surface, then 
a flow expansion occurs at the turning angle region, until the 
maximum local Mach number is reached. A weak shock is formed 
at the end of the boattail region as shown in Figure (24), 
while for the viscous flow, the weak shock develops away from 
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boattail surface due to the boundary layer and the separation 
of the flow. The circulating flow in the base region is -more 
significant than in the inviscid flow as depicted in Figure 
(25) . In addition, from Table 2, for viscous flow the maximiam 
local Mach number for the conical afterbody is higher than for 
the contoured afterbody. 

For each afterbody the pressure drag coefficient (C^j) 
versus free stream Mach number (M^^) for inviscid and viscous 
flow are plotted in Figure (26) . The negative sign of is 
due to the fact that the calculation of pressure starts from 
mid-section through the base of afterbody and ignores the 
forebody pressure. These results show that the pressure drag 
coefficient is higher for viscous flows than for inviscid 
flows for each given afterbody. This may be caused by the 
viscosity effect and the pressure distribution difference in 
the flow field. In addition, it can be seen from the chart, 
that the for the afterbody with small perturbation solution 
contour Mach l.lO (SPSl.l) has the lowest values over the 
entire Mach number range. Therefore, the SPSl.l contour shows 
to be relatively the best among these afterbodies. 

b. Complete Bodies 

Complete bodies consist of an arbitrary forebody 
(ellipsoid) joined to various afterbodies such as the small 
perturbation solution contours (SPS_1.1 and SPS_1.2) and a 
conical afterbody. The approaching free stream Mach number (M„ 
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ranging from 0.50 to 1.50) starts from outer boundary with 
distance of 1 times for the coarse grid and 5 times body's 
length for the fine grid respectively. Significant differences 
between coarse and fine grid calculations can be seen from the 
characteristic bow shock. In the coarse grid, the bow shock 
hits the outer boundary, this causes the approaching free 
stream Mach number not to be the same as in the input 
parameter as shown in Figure (27) . While in the fine grid, the 
bow shock dies out before reaches the outer boundary as 
depicted in Figure (28) . So, further discussion is focused on 
the fine grid exclusively. 

The flow stagnates on the nose and then follows the body 
surface until it reaches a maximum local Mach number as 
tabulated in Table 3 . The critical Mach number for these 
complete bodies is approximately at M„ = 0.70, where the 
maximum local Mach number reaches unity at the shoulder 
region. 

Figure (29) shows a typical high subsonic free stream Mach 
number (M„, = 0.85) flow over the complete body SPS_1.1. The 
flow stagnates on the nose tip, then flow is accelerated 
following the forebody surface reaching a sonic line at the 
mid- way of the forebody and forming a supersonic region at the 
shoulder. Then, the flow deccelerates at the mid- section and 
accelerates again until reaches a maximum local Mach number at 
the starting region of the afterbody and a weak shock occurs 
in this region. A supersonic approaching free stream Mach 
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number, is depicted in Figures (30) and (31) for SPS_1.1 and 
SPS_1.2 respectively. A bow shoclc is obviously seen in front 
of the nose. The characteristic of the bow shoclc is more 
inclined down stream for the higher as shown in Figure 
(32) . A subsonic region is formed between bow shoclc and the 
nose, then the flow accelerates along the forebody surface up 
to a supersonic region in the mid- section. The expansion flow 
occurs in the starting contour region until it reaches a 
maximum local Mach number. Then, a shock is formed in the 
contour region. Similar as in the afterbody only, the shock 
location depends upon the contour surface and M„. The shock 
location for a given contour is more downstream for higher 
and at the same M^g, the shock location for SPS_l.l is more 
downstream than SPS_1.2. A typical residual history of CFD 
calculation for complete body with small perturbation solution 
contour is shown in Figure (33) . Convergence is obtained at 
about 500 iterations. 

For the conical afterbody, the flow characteristic is the 
same as the other complete bodies up to the mid- section 
region. The expansion flow occurs at the turning angle, then 
the flow accelerated along the conical surface and weak shock 
is formed at the edge of base. Figures (34) and (35) show the 
Mach number contour and corresponding residual calculation for 
conical afterbody at M„g = 1.10. 

The pressure drag coefficient (C^j-press) versus free 
stream Mach number (M^,) for fine grid complete bodies are 
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plotced in Figure (36) . The drag rises sharply in the high 
subsonic Mach number (M„ = 0.95) and reaches a maximum (peak) 
at M„ = 1.10. Then the drag decreases with a shallow curve as 
the M„ increases. The decreasing shallow curve may be caused 
by the bluntness of the nose and it agrees with Shapiro 
[Ref. 18] because the fineness ratio and bluntness of the nose 
of bodies of revolution are the important factors that 
contribute the drag curve at transonic and supersonic range. 
It can be seen from the graph, the drag curve of complete body 
with conical afterbody is higher than with small perturbation 
solution contour at the entire M„. Furthermore, the peak of 
the drag curve is approximately 15% higher. The drag curve for 
SPS_1.1 and SPS_1.2 are likely to have the same trend up to M„ 
= 0.95; beyond this Mach number, the drag curve for SPS_1.2 is 
slightly greater than for SPS_1.1. Therefore, the complete 
body with small perturbation solution contour afterbody Mach 
1.1 (SPS_1.1) relatively gives the lowest drag. 
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VI. CONCLUSIONS AND RECOMMENDATIONS 



The use of numerical simulation (CFD) appears to be the 
most cost effective method to predict the aerodynamic 
performance, especially in the transonic range. 

In this research, the grid-generating program GRAPE is 
suitable only for the geometry with smooth contour surface 
(SPS_1.1 and SPS_1.2) , while the program GRIDGEN2D is used for 
the geometry of a conical afterbody with non- smooth points. 

We have shown that for a complete body model the use of a 
fine grid (77x9x120) is more reliable than a coarse grid 
(77x13x60) . This is shown by the characteristic of the bow 
shock at > 1 as seen in Figures 27 and 28. The pressure 
drag coefficient (C^j) versus free stream Mach number (M„) 
graphs show that the small perturbation solution contour for 
Mach 1.10 (SPS_1.1) gives relatively the lowest on both 
models (afterbody and complete body) , a decrease of 
approximately 15% of peak in the transonic range compared to 
the conical afterbody. Therefore, the best design of an axi- 
symmetric body such as missiles, projectiles and aircraft, can 
be based on the small perturbation solution contour. 

Finally, this work may be continued with the investigation 
a complete body with pointed nose geometry and the application 
of a Euler as well as Navier-Stokes flow-solvers within 
various small angle of attack. 
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TABLE 1. INPUT PARAMETER FOR OVERFLOW. IN 
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inflow/ outflow. 



TABLE 2. MAXIMUM LOCAL MACH # FOR AFTERBODIES 
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TABLE 3. MAXIMUM LOCAL MACH NUMBER FOR COMPLETE BODY 
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KSI(30 vs. X For M.= 1.05,1.10,1.20 
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Figure 2. Numerical solution of equation (33). 
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Figure 4. Numerical integration of equation (63) (normalized) 
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Rgvire 6. 2D grid , afterbo^ model with SPS_1,1 
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Rgure 7. 2D grid, complete bo<fy with SPS_L1 . 
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Figure 8. 3D grid , afterbo^ model with SPS l.l 
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Rgure 9. 3D grid , complete-body with SPS1.1 
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Figure 10. 2D grid, conical afterbody model. 
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Figure 11. 2D grid, complete body with conical afterbody. 
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Figure 12. 3D grid, conical afterbody model 
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Figure 13. 3D grid , complete body with conical afterbody . 
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Figure 14. 3D fine-grid, complete body with conical afterbody. 
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Figure 15. 3D fine-grid, complete body with SPS_1.1 . 
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Figure 16. Afterbody model SPS l.l , invisdd flow , Macb 1.10 . 
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Figure 17. Afterbody model SPS l.l , iuvisdd flow , Macfa 1.20 . 
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Figure 18. Afterbocfy model SPS1.2 , invisdd flow , Mach 1.10 , 
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Rgure 19. Afterbody model SPS_1.2 , inviscid flow , Mach 1.20 . 
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Figure 20. Afterbody model SPS l-0 , inviscid flow , Mach 1.10 . 
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Figure 22. Afterborfy model SPS l.l , viscous flew , Mach 1.10 . 
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Figure 23. Afterbody model SPS1.2 , viscouc flow , Mach 1.20 . 
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Rgvire 24. Conical afterbody model , inviscid flow , Mach 1.10 . 
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figure tJorncaTatterbody model , viscous flow , Madi 1.10 - 



Cd-press vs Mach # for afterbody models 
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Figure 26. ^ - press, versus Mach number for afterbody models. 
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Figure 27. Coarse-grid (77x13x60) for half complete body models. 
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Figure 29. Quarter complete body aft-SPS_l.l , Mach 0.85 . 
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Figure 30. Quarter complete body aft-SPS_l.l , Mach 1.10 . 
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Figure 31. Quarter complete body aft-SPS_1.2 , Mach 1.20 . 
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Figure 32. Quarter complete body aft-SPS_l.l , Mach 1.20 . 
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Figure 33. Residual history of complete body SPS_1.1 , Mach 1.10 . 
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Figure 34. Quarter complete body aft-conical , Mach 1.10 . 
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Figure 35. Residual history of complete bcxly aft-conial , Mach I.IU . 
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Figure 36. C . - press, versus Mach number for quarter complete body models . 
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PROGRAM KSI.F 

************************************************************* 

* This program is developed to calculate KSI (X) , K(X) and X * 

* using mumerical integration based on trapezoidal rule to * 

* solve equation (32) and plotting the output as shown in * 

* Figs . 1 . * 



REAL M(3) , X(0;401,3), A, A1 , B(401,3) 

# , KSI, P, FUNC, H, DD, Q, L, K(401,3), N, N1 

INTEGER YY 

OPEN (UNIT=9 , FILE=' KSI' , STATUS =' UNKNOWN' ) 

PRINT *, 'ENTER LOWER & UPPER BOUND, # OF INTERVALS, # OF 
#DATA SET' 

READ *, A, Al, N, N1 
100 DO 10 1=1,3 

PRINT* , ' ENTER MACH NO . ' 

READ*, M(I) 

C M(I) = 0. 9+1*0. 1 

IF (M(I) .LT.1.0) THEN 
P = -1. 

ELSE 
P = 1. 

ENDIF 

DO 20 J=1,N1 

B(J,I) = J*A1/N1 
H = (B(J,I) -A) /N 
AREA = 0. 

K(J,I) = (-0.0102/(M(I) **2) )*( (ABS( (B(J,I)**2) - ( (ABS (1- 
# M(I) **2) ) **1.7574) ) ) ** (2./3 . ) + (ABS (l-M(I) **2) ) **1.1716) 
C 

DO 30 L= 1,N 

KSI = A + (L-0.5) *H 

DD = (KSI**2+P* (ABS (l-M(I) **2) ) **1.7574) 

IF (DD.GT.0.0) THEN 
FUNC = 1./ (DD** (1./3 . ) ) 

Q = 1. 

ELSE 

FUNC = 1. / ( ( -1. *DD) ** (1. /3 . ) ) 

Q = -1. 

ENDIF 

AREA = AREA + H*FUNC*Q 
X(J,I) = AREA 
30 CONTINUE 

20 CONTINUE 

10 CONTINUE 

C 

DO 40 J= 1,N1 

PRINT 50, (B(J,I) ,X(J,I) ,K(J,I) ,1=1,3) 
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50 F0RMAT(1X,3 (F6.3,2F10.5,1X) ) 

40 CONTINUE 

200 PRINT*, 'TRY AGAIN?, ENTER 1 FOR YES , 2 FOR DATA FILE, 
# OTHERS FOR NO' 

READ*, YY 
IF (YY.EQ.l) THEN 
GO TO 100 
ELSE 

IF (YY.EQ.2) THEN 
DO 110 J= 1,N1 

WRITE (9,50) (B(J,I) ,X(J, I) ,K(J, I) ,1=1,3) 

110 CONTINUE 

GO TO 200 
ENDIF 

ENDIF 

END 



PROGRAM ZETA.F 

************************************************************ 

* This program is written to calculate ZETA(r) and r using* 

* equation (29) and plotting output in Fig. 3 . * 

★★★★★★*★★★★★★******★**★★**★★★*★★★★★★★★★★***★** 

REALM(3), ZETA(0:14,3) , R(14) 

OPEN(UNIT=9, FILE=' ZZ' , STA'TUS= ' UNKNOWN' ) 

C 

DO 10 1= 1,3 

C M(I) = 0. 7+1*0. 1 

PRINT* , ' ENTER MACH NO . ' 

READ*, M(I) 

DO 20 J= 1,14 
ZETA(0,I) = 1000.0 
PRINT*, 'ENTER R VALUES' 

READ*, R(J) 

C R(J) = 0.1*J 

ZETA(J,I) = (4/ (R(J) **2) ) + (ABS (1- (M(I) **2) ) 

# * (R(J) **2.8284) ) 

# +( (l-M(I) **2) *(R(J)**7.657) /50.63 

IF (ZETA(J, I) .GE.ZETA(J-1, I) ZETA ( J, I) =ZETA ( J- 1 , I) 

20 CONTINUE 

10 CONTINUE 
C 

PRINT 30, M 

30 FORMAT(11X,5F10.4) 

DO 40 J= 1,14 

WRITE(9,50) (R(J), ( ZETA ( J , I ) , 1=1 , 4 ) 

50 FORMAT(4X,3 (2F10.4,2X) 

40 CONTINUE 
END 
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PROGRAM D2D3.F 

*********************************************************** 

* This program will read data files output Gridgen2d * 

* from 2D to 3D for further processing with rotategr.f * 

REAL X(77,60) ,Y(77,60) ,Z(77,60) 

READ (12,*) IDIM,JDIM 

READ ( 12 , * ) ( (X ( I , J) , 1=1 , IDIM) , J=1 , JDIM) , 

# ( (Y(I, J) ,I=1,IDIM) , J=l, JDIM) 

Z(I, J)=0.0 

WRITE (14,*) IDIM, JDIM, 1 

WRITE (14,*) ( (X(I, J) ,1=1, IDIM) , J=1,JDIM) , 

# ( (0.0, 1=1, IDIM) ,J=1, JDIM) , 

# ( (Y(I, J) ,I=1,IDIM) , J=l, JDIM) 

STOP 

END 



PROGRAM HALF.F 

* This program will write out half of a grid file output * 

* from Grape for further processing with rotategr.f * 



DIMENSION X(200,200) , Y(200,200) 

READ (30,*) IDIM, JDIM 

READ (30,*) ( (X(I,J) ,1=1, IDIM) , J=l, JDIM) , 

# ( (Y(I, J) , I=1,IDIM) , J=l, JDIM) 



c 

IDMD2=IDIM/2 + 1 
C 

DO 10 J=1,JDIM 

Y(IDIM2,J) =0.0 
10 CONTINUE 



C 



WRITE (31, *) 
WRITE (31, *) 

# 

# 

STOP 

END 



IDMD2, JDIM, 1 

( (X(I, J) , I=IDMD2) , J=l, JDIM) , 
((0.0, I=IDMD2 ) , J=1 , JDIM) , 

( (Y(I, J) , I=IDMD2) , J=l, JDIM) 
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PROGRAM ROTATEGR.F 

**★★★★★★★**★*★★★*★★★*★*★★*★*★*★**★★★★*★★***★★★★*****★★★★★★★★★ 

* This program will create 3D grid (body of revolution) * 

* by rotating 3D plane grid for further processing with * 

* Flow Solver * 

************************************************************* 



c 
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DIMENSION X (200, 200) , Y(200,200), Z(200,200) 
DIMENSION XX(200, 100,200) , YY (200 , 100 , 200) , 
#ZZ (200,100,200) 

CHARACTER* 30 FNI 
CHARACTER* 30 FNO 



PRINT*, 'INPUT GRID' 

READ (*,21) FNI 

PRINT*, 'SG =0, MG = 1', ',UNF=1, FORM=2 ' 

READ ( * , * ) IGRI 

PRINT*, ' IFORI = 1 UNF, IFORI = 2 FORM' 
READ ( * , * ) IFORI 
FORMAT (A) 

IF ( IFORI -EQ.l) THEN 
REWIND 1 

OPEN (1, FILE= FNI, FORM=' UNFORMATTED ' ) 

IF (IGRI. EQ.l) READ (1) MGR 
READ (1) I1,J1,K1 

READ (1) ( (X(I,J) ,1=1,11) ,J=1,J1) , 

# ( (Y(I,J) ,1 = 1,11) ,J=1,J1) , 

# ( (Z(I,J) ,1=1,11) ,J=1,J1) 



CLOSE (1) 

ENDIF 

IF (IFORI.EQ.2) THEN 
REWIND 2 

OPEN (2 , FILE=FNI, FORM=' FORMATTED' ) 
IF (IGRI. EQ.l) READ (2,*) MGR 
READ (2,*) I1,J1,K1 

READ (2,*) ( (X(I, J) ,1=1,11) , J=l, Jl) , 



# ( (Y(I,J) ,1=1,11) ,J=1,J1) , 

# ( (Z(I,J) ,1=1,11) ,J=1,J1) 



CLOSE (2) 
ENDIF 



PI = 4.*ATAN(1.) 

PRINT*, 'NO OF PLANE IN J DIRECTION ?' 
READ ( 5 , * ) JM 

DTH = (180./ (JM-1) )* (PI/180) 

DO 11 1=1,11 
DO 11 J=1,J1 
K=J 

XX(I,2,K) = X(I,K) 

YY(I,2,K) = 0.0 
ZZ(I,2,K) = Z(I,K) 

11 CONTINUE 
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c 



20 



c 



30 



IM=I1 

KM=J1 

DO 20 J=3,JM+1 
DO 20 1=1, IM 
DO 20 K=1,KM 
XX(I, J,K) =X(I,K) 

TH = (J-1)*DTH 
YY(I,J,K) = SIN(TH) *Z (I, K) 
ZZ(I,J,K) = COS (TH) *Z(I, K) 
CONTINUE 



DO 30 1=1, IM 
DO 30 K=1,KM 



J=1 

XX(I, J,K) 
YY(I, J,K) 
ZZ (I, J,K) 
J=JM+2 
XX(I, J,K) 
YY(I, J,K) 
ZZ (I, J, K) 
CONTINUE 



= XX(I,J+2,K) 
=-YY(I, J+2,K) 
= ZZ(I,J+2,K) 

= XX(I,J-2,K) 
=-YY(I, J-2,K) 
= ZZ(I,J-2,K) 



JM=JM+2 



PRINT*, 'Output filename =' 

READ (5,21) FNO 
C 

REWIND 3 

OPEN( 3, FILE=FNO, FORM= ' UNFORMATTED ' ) 

WRITE (3) IM,JM,KM 

WRITE(3) (((X(I,J,K), 1=1, IM) , J=l, JM) ,K=1,KM) , 

# (((Y(I,J,K), 1=1, IM) , J=l, JM) ,K=1,KM) , 

# (((Z(I,J,K), 1=1, IM) , J=l, JM) ,K=1,KM) 
CLOSE (3) 

STOP 

END 
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PROGRAM READX.F 

************************************************************ 

* This program will read 3D grid file and convert from * 

* formatted into unformatted file (grid. for to grid. in) * 

* for further processing with Overflow. * 

************************************************************ 

DIMENSION X( 77, 13,60) , Y(77,13,60), Z(77,13,60) 

OPEN (UNIT=12, FILE=' GRID. FOR' , STATUS =' UNKNOWN' ) 

OPEN (UNIT=14,FILE=' GRID. IN' , STATUS= ' NEW' , 

#FORM= ' UNFORMATTED ' ) 

READ (12,*) IDIM, JDIM,KDIM 

READ(12,20) (((X(I,J,K), 1=1 , IDIM) , J=1 , JDIM) , K=1 , KDIM) , 

# (((Y(I,J,K), 1=1, IDIM) ,J=1, JDIM) ,K=1, KDIM) , 

# (((Z(I,J,K), 1=1, IDIM) ,J=1, JDIM) ,K=1, KDIM) 
c 

WRITE (14) IDIM, JDIM, KDIM 

WRITE(14) (((X(I,J,K), 1=1, IDIM) ,J=1, JDIM) ,K=1, KDIM) , 

# (((Y(I,J,K), 1=1, IDIM) ,J=1, JDIM) ,K=1, KDIM) , 

# (((Z(I,J,K), 1=1, IDIM) ,J=1, JDIM) ,K=1, KDIM) 

20 FORMAT (6F15. 6) 

STOP 

END 



PROGRAM READQ.F 

************************************************************ 

* This program will read q.save file output from * 

* Overflow and convert from unformatted into formatted * 

* file for further processing with Plot3d. * 

************************************************************ 



DIMENSION 0(77,13,60,5) 

OPEN (UNIT=1, FILE=' Q.SAVE' , STATUS= ' OLD ' , 

#FORM= ' UNFORMATTED ' ) 

OPEN(20, FILE='Q.FORM' , STATUS= ' NEW' , F0RM= ' FORMATTED ' ) 

READ(l) NI,NJ,NK 

READ(l) FSMACH,ALPHA,RE,TIME 

READ(l) ( ( ( (0(1, J,K,NX) ,I=1,NI) , J=1,NJ) ,K=1,NK) ,NX=1,5) 

WRITE (20,*) NI,NJ,NK 

WRITE (20,*) FSMACH, ALPHA, RE, TIME 

WRITE(20,*) ( ( ( (0(1, J,K,NX) 



#,I=1,NI) , J=1,NJ) ,K=1,NK) ,NX=1,5) 
PRINT*, NI,NJ,NK 
STOP 
END 
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PROGRAM OVERFLOW. IN 

icicicieie'kieieicicicicieieicicieieieiciricicicicicieicicicicicic'k'kic'kic'kie'k'k'k'k'k'k'k'k'k'k'k'k'k'k'kie'k'k'k'k 

* Input file for Overflow to run either EULER or NAVIER- * 

* STOKES Flow- solver. * 

$GLOBAL 

CHIMRA=.F., NSTEP=100, RESTART=.F., NSAVE=100, NQT=0, 
$END 
$FLOINP 

ALPHA=180.0, FSMACH=0.95, REY= 6.00E6, TINF= 520.0, 
$END 
$GRDNAM 

NAME= 'AXI- SYMMETRIC BODY WITH AFT- CONE, 77x13x60 
GRID' , 

$END 

$METPRM 

IRHS = 0, ILHS = 2, IDISS = 2, 

$END 

$TIMACU 

DT = 0.1, ITIME = 1, TFOSO = 1.00, CFLMIN = 5.00, 

$END 

$SMOACU 

ISPECJ = 2, DIS2J = 2.00, DIS4J = 0.02, 

ISPECK = 2, DIS2K = 2.00, DIS4K = 0.02, 

ISPECL = 2, DIS2L = 2.00, DIS4L = 0.02, 

SMOO = 1.00, EPSE = 0.35, 

$END 

$VISINP 

VISCJ = .F. , VISCK = .F. , VISCL = .F. , 

NTURB = 0, 

ITTYP = 1 , 

ITDIR = 3 , 

JTLS = 1 , 

JTLE = 77, 

KTLS = 1 , 

KTLE = 13 , 

LTLS = 1 , 

LTLE = 60, 

TLPAR1= .3, 

$END 

$BCINP 

NBC = 6, 



IBTYP 


= 


15, 


12, 


12, 


1, 


32, 


15, 


IBDIR 


= 


1, 


2, 


-2, 


3, 


3, 


-1, 


JBCS 


= 


1, 


1, 


1, 


1, 


1, 


77, 


JBCE 


= 


1, 


77, 


77, 


77, 


77, 


77, 


KBCS 


=: 


1, 


1, 


13, 


1 , 


1, 


1, 


KBCE 


= 


13, 


1, 


13, 


13, 


13, 


13, 


LBCS 


= 


1 , 


1 , 


1, 


1, 


60, 


1, 


LBCE 


= 


60, 


60, 


60, 


1, 


60, 


60, 



$END 
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